%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file computes the initial steady state of the baseline model

%% Legend
% IPF: Adrian, Erceg, Kolasa, Linde, and Zabczyk (IMFwp 2021) (Integrated Policy Framework)
% DICE: Barrage and Nordhaus (2023)
% CHM: Carattini, Heutel, and Melkadze (RED forth)
% HKS: Hassler, Krusell, Smith (Macro Hand. 2016)
% HKO: Hassler, Krusell, Olovsson (PIE 2024)
% EHMS: Ernst, Hinterlang, Mahle, Stähler (JIE 2023)
% LM: Liberati and Marinelli (QEF 2021)
% WEO: World Economic Outlook, october 2023
% OWD: Our world in data
% H: Home (advanced economy)
% F: Foreign (emerging economy)

clear all; 
clc; 
close all;

SIMUL=1;   % if 0 you want to run console only, if 1 you run also console_final_opt
SYM=0;     % if 1 calibration is symmetric, if 0 the calibration is asymmetric


%% Parameters and targets constant across calibrations
% Parameters
psiX=0.5;             % pollution convexity
eme_factor_dis=1;   % how much EMEs suffer from pollution wrt to AE (eme_factor_dis=2, HKO)
eme_factor_ab=1;    % how much abatment is costly in EME wrt to AE (eme_factor_ab=1.89, EHMS)
delta=0.1;          % depreciation rate (standard)
ome=0.3;            % trade openness (IPF)
eta=1.5;            % elasticity of substitution between domestic and foreign goods (IPF)
iota=1.0135;        % labor productivity growth (to have a good match with DICE)
chi=1;              % inverse of Frisch elasticity
nu=1.6;             % convexity of the abatement cost (CHM)
csi=2;              % elasticity of substitution between green and brown goods (CHM)
zeta=0.332;         % share of brown output (CHM)
GAM=0.02;           % H shallowness of financial markets (IPF) 
deltax=1-(0.9965^4);% pollution depereciation in the Heutel framework (CHM)
phi=0.0023;         % decay rate of emissions that depreciate, but not within a decade (this share is (1-phiL)*phi0) (HKS)
phi0=0.401;         % (1-phi0) is the share of depreciating emissions that fade away in one year (HKS)
phiL=0.2;           % share of emissions that do not depreciate (HKS)
Asym=1;             % H TFP in the symmetric calibration (noralization)
mu=0.05;            % H abatement rate (DICE)
muz=0.05;           % F abatement rate (DICE)
price_clean=514;    % carbon price with zero emission (DICE)

% Targets
r=1.025;            % H real interest rate (IPF)
gY=0.14;            % H public expenditure GDP ratio  (IPF)
gYz=0.14;           % F public expenditure GDP ratio (IPF)
iY=0.22;            % H investment GDP ratio (World Bank - World Development Indicators)
FDIY=0.086;         % H stock of FDI outflows GDP ratio (IMF CDIS)
xGtC=887;           % 2020 level of atmospheric carbon in GtC  (DICE) 
xbar=581;           % pre-industrial level of carbon in GtC (HKS)
x1_base=684;        % 2020 permanent carbon sink (HKS)
eWGtC=9.54;         % 2020 global emissions in GtC          (OWD)
gdp_world=87325;    % world GDP in dollar billion 2019      (WEO)

%% Parameters and targets different across calibrations
if SYM==1
% Parameters    
n=0.5;              % H population share (symmetric)
lam=1;              % H-F GDP per capita ratio (symmetric) 
GAMz=GAM;           % F shallowness of financial markets (symmetric) 

% Targets  
rz=r;               % F real interest rate EME (symmetric)
bY=0;               % H NFA position (symmetric)
iYz=iY;             % F investment GDP ratio (symmetric)
FDIYz=FDIY;         % F stock of FDI outflows GDP ratio (symmetric)
eGtC_share=0.5;     % H share of emissions  (symmetric)

else
% Parameters
n=0.14;             % H population share  (2020 value, WEO)  
lam=4.47;           % H-F GDP per capita ratio (2020 value, WEO)                   
GAMz=0.06;          % F shallowness of financial markets (IPF) 

% Targets
rz=1.03;            % F real interest rate (IPF)
bY=0.22;            % A NFA position (IPF)
iYz=0.33;           % F investment GDP ratio (World Bank - World Development Indicators)
FDIYz=0.045;        % F stock of FDI outflows GDP ratio (IMF CDIS)
eGtC_share=0.294;   % H share of emissions 0.294      (2020 value, OWD)


end


%% Ex-post parameters and steady-state values
gamma= ome*(1-n)/(lam*n+1-n);           % AE import quasi-share
gammaz= ome*n*lam/(n*lam+1-n);          % EME import quasi-share
beta=iota/r;                            % AE discount factor
betaz=iota/rz;                          % EME discount factor
rkB=iota/beta-(1-delta);            % brown return of AE domestic capital
rkB_zd=iota/beta-(1-delta);         % brown return of EME foreign capital
rkB_zu=iota/betaz-(1-delta);        % brown return of AE foreign capital
rkB_zz=iota/betaz-(1-delta);        % brown return of EME domestic capital
rkG=iota/beta-(1-delta);            % green return of AE domestic capital
rkG_zz=iota/betaz-(1-delta);        % green return of EME domestic capital
rkG_zd=iota/beta-(1-delta);         % green return of EME foreign capital
rkG_zu=iota/betaz-(1-delta);        % green return of AE foreign capital



options = optimoptions('fsolve','MaxFunEvals',300000,'MaxIter',30000,'TolFun',1e-15);
load gdpW_sym.mat
x0=[1.1698    1.0000    1.0044    1.0044 27.5878 27.5878 1];
x0=x0+0.001;

X = fsolve(@(XX) find_steady(XX,gamma,gammaz,eta,lam,bY,beta,betaz,GAM,GAMz,n,gY,gYz,FDIY,FDIYz,delta,iota,...
                                iY,iYz,zeta,csi,mu,muz,nu,rkB,rkG,rkB_zu,rkB_zd,rkG_zu,rkG_zd,rkB_zz,chi,...
                                rkG_zz,r,rz,price_clean,gdp_world,eWGtC,eGtC_share,Asym,gdp_sym,SYM,eme_factor_ab),x0,options);
yH=X(1);                                              
pH=X(2);
pB=X(3);
pBz=X(4);
zetae=X(5);
zetaez=X(6);
A=X(7);


x0=real(X);
save x0 x0;

pF=((1-(1-gamma)*pH^(1-eta))/gamma)^(1/(1-eta));
s=(gammaz*pH^(1-eta)+(1-gammaz)*pF^(1-eta))^(1/(1-eta));
gdp=pH*yH;
yFz=gdp/(lam*pF);
gdpz=pF*yFz/s;
b=bY*gdp;
dH=(1-beta/betaz)/GAM;
dF=s/GAMz*(betaz/beta-1);
v=-(b+dH+dF*(1-n)/n);
g=gY*yH;
gz=gYz*yFz;
k_zd=n*FDIY*gdp/(s*(1-n));
k_zu=FDIYz*gdpz*s*(1-n)/n;
i_zd=k_zd*(1-(1-delta)/iota);
i_zu=k_zu*(1-(1-delta)/iota);
i=iY*gdp-i_zu;
i_zz=iYz*gdpz-i_zd;
k=i/(1-(1-delta)/iota);
k_zz=i_zz/(1-(1-delta)/iota);
yB=zeta*(pB/pH)^(-csi)*yH;
conv_gdp=gdp_world/(n*gdp+(1-n)*pF*yFz);  
kappaA=price_clean*3.67/(conv_gdp*(n/zetae+(1-n)*s*eme_factor_ab/zetaez));
kappaAz=kappaA*eme_factor_ab;
tau=(mu^nu)*kappaA/zetae;                    
tauz=(muz^nu)*kappaAz/zetaez;                  
pB2=pB-tau*(1-mu)*zetae-kappaA/(1+nu)*mu^(1+nu);
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);
pGY=pH*yH-pB*yB;
alpha1=k/(iota*(pB2/rkB*yB+pGY/rkG));
kB=alpha1*iota/rkB*pB2*yB;
kG=alpha1*iota/rkG*pGY;
alpha2=k_zu/(iota*(pB2/rkB_zu*yB+pGY/rkG_zu));
kB_zu=alpha2*iota/rkB_zu*pB2*yB;
kG_zu=alpha2*iota/rkG_zu*pGY;
hB=(yB/(A*(kB/iota)^alpha1*(kB_zu/iota)^alpha2))^(1/(1-alpha1-alpha2));
w=(1-alpha1-alpha2)*pB2*yB/hB;
hG=(1-alpha1-alpha2)*pGY/w;
h=hG+hB;
yG=A*(kG/iota)^alpha1*(kG_zu/iota)^alpha2*(hG)^(1-alpha1-alpha2);
pG=pGY/yG;
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;
pGYz=gdpz-pBz*yBz;
alpha1z=k_zz/(iota*(pB2z/rkB_zz*yBz+pGYz/rkG_zz));
kB_zz=alpha1z*iota*pB2z*yBz/rkB_zz;
kG_zz=alpha1z*iota*pGYz/rkG_zz;
alpha2z=k_zd/(iota*(pB2z/rkB_zd*yBz+pGYz/rkG_zd));
kB_zd=alpha2z*iota*pB2z*yBz/rkB_zd;
kG_zd=alpha2z*iota*pGYz/rkG_zd;
c=-(i+pH*g+kappaA/(1+nu)*mu^(1+nu)*yB+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd    +rkG_zd*kG_zd)    /iota+(rkB_zu*kB_zu    +rkG_zu*kG_zu)    /iota+rz            *(dH+v)    /iota-r/iota    *(b+    dH    +v)-gdp);
cz=n/((1-n)*gammaz)*(pH/s)^(eta)*(yH-((1-gamma)*pH^(-eta)*(c+i+kappaA/(1+nu)*mu^(1+nu)*yB+i_zu)+g))-(i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
wz=(cz^(1/chi)*(1-alpha1z-alpha2z)*(pB2z*yBz+pGYz))^(chi/(1+chi));
hz=(wz/cz)^(1/chi);
hBz=(1-alpha1z-alpha2z)*pB2z*yBz/wz;
hGz=hz-hBz;
Az=yBz/((kB_zz/iota)^alpha1z*(kB_zd/iota)^alpha2z*(hBz)^(1-alpha1z-alpha2z));
yGz=Az*(kG_zz/iota)^alpha1z*(kG_zd/iota)^alpha2z*(hGz)^(1-alpha1z-alpha2z);
pGz=pGYz/yGz;
e=zetae*(1-mu)*yB;                                                                                                 
ez=zetaez*(1-muz)*yBz;  

bz=-n*b/((1-n)*s);
vz=-v/s;
dHz=-dH/s;
dFZ=-dF/s;
mpz=pH/s*gammaz*(pH/s)^(-eta)*(cz+i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
xpz=n/(1-n)*pF/s*gamma *pF ^(-eta)*(c+i +i_zu+kappaA/(1+nu)*mu^(1+nu)*yB);
tbz=xpz-mpz;
prova=cz+i_zu+i_zz+pF/s*gz+kappaAz/(1+nu)*muz^(1+nu)*yBz+tbz-(pF/s*yFz);

tau_clean=kappaA/zetae;
tauz_clean=kappaAz/zetaez;

%% Other parameters
kappaX=(5.3*10^(-5)*(1+psiX)*(xGtC-xbar)^(-psiX))/(n+(1-n)*eme_factor_dis);         % pollution distulity (as in GHKT)

kappaXz=eme_factor_dis*kappaX;

%% Define steady state values
yHss=yH;
yFzss=yFz;
pBss=pB;
pBzss=pBz;
pHss=pH;
bss=b;
kBss=kB;
k_zdss=k_zd;
k_zuss=k_zu;
rkGss=rkG;
rkG_zzss=rkG_zz;
tauss=tau;
tauzss=tauz;
iss=i;
i_zzss=i_zz;
css=c;
czss=cz;
muss=mu;
muzss=muz;
zetaess=zetae;

if SYM==1
save gdpW_sym gdp_sym
end


%% Save parameters
save par gamma gammaz eta beta betaz GAM GAMz n delta deltax iota zeta csi kappaA kappaAz nu v g gz alpha1 alpha2 alpha1z...
     alpha2z A Az zetaess conv_gdp zetaez x1_base yHss yFzss pBss pBzss pHss bss kBss k_zdss k_zuss rkGss chi SYM...
     rkG_zzss tauss tauzss iss i_zzss css czss muss muzss kappaX psiX  xbar tau_clean tauz_clean phi phi0 phiL kappaXz

if SIMUL==1
console_final_opt
end